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strong to weak coupling behaviour. Strong correlations are evident in the weak-coupling 
regime. 
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1 Introduction 

While there has been considerable recent interest in graphene sparked by its discovery 
and subsequent experimental study [TJ, the remarkable properties of electronic systems 
on two-dimensional honeycomb lattices have been suspected for many years [2]. In brief, 
for a carbon monolayer having one mobile electron per atom, a simple tight-binding 
model shows that the spectrum of low-energy excitations exhibits a linear dispersion 
relation centred on zeroes located at the six corners of the first Brillouin zone (eg. 
[3]). Using a linear transformation among the fields at two independent zeroes it is 
possible to recast the Hamiltonian in Dirac form with Nf = 2 flavors of four-component 
spinor ip, the counting of the massless degrees of freedom coming from 2 C atoms per 
unit cell x 2 zeroes per zone x 2 physical spin components per electron. Electron 
propagation in the graphene layer is thus relativistic, albeit at a speed vp ~ c/100. 
The implications for the high mobility of the resulting charge carriers (which may be 
negatively-charged "particles" or positively-charged "holes" depending on doping) is the 
source of the current excitement. The stability of the zero-energy points is topological 
in origin, as emphasised in [I]. 
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While the above considerations apply quite generally, a realistic model of graphene 
must incorporate interactions between charge carriers. One such model due to Son 
[5] has Nf massless fermion flavors propagating in the plane, but interacting via an 
instantaneous 3d Coulomb interaction. In Euclidean metric and static gauge doAo = 
the action reads 



N f \ r 

Si = y^l dx d 2 x(ifj a 'yodo^a + v F ifj a j.Vif} a + iV^al^a) + ^2 / dx d s x(diVf, (1) 

a=l 



where e is the electron charge, V = Aq is the electrostatic potential, and the 4x4 Dirac 
matrices satisfy {7^,7*,} = 25^ u , \i = 0, 1,2,3. In our notation x is a vector in the 2c? 
plane while the index % runs over all three spatial directions. Within the graphene layer, 
classical propagation of the potential is obtained by integrating over the perpendicular 
coordinate, yielding 



e- 



DM = m (2) 

To proceed, assume a large-iV/ limit so that the dominant quantum correction Tl(p) 
comes from a vacuum polarisation fermion - antifermion loop. The resummed V prop- 
agator becomes 

BlW = (^ W -n W) -=(f + ^)"', (3) 

where p 2 = (po,p) 2 = pi + v 2 F \p\ 2 . In either the strong coupling or large- Nf limits Di(p) 
is thus dominated by the quantum correction, the relative importance of the original 
Coulomb interaction being governed by a parameter A = |n/_D | po=0 . Restoring SI 
units, we obtain 

e 2 N f 

A = 7 ~ 1.7N f . (4) 
16e hv F } y ' 

The form of the interaction ([3]) means that analytic methods are trustworthy in 
the large-iVj regime. For instance, in the strong coupling limit e — ► we expect 
a modification of the dispersion relation, such that the fermion energy is related to 
momentum via uj oc p z , where z is a dynamical critical exponent predicted to take the 
value z ~ 1 — ^57^ ~ 0.8 for Nf = 2 [5]. Ref. [S] in addition discusses the phase 

diagram of the graphene model (JTJ) in the (Nf, e~ 2 ) plane, and raises the possibility 
of symmetry breaking due to non-perturbative Nj 1 effects. The symmetry breaking, 
due to the spontaneous condensation of particle - hole pairs, is signalled by an order 
parameter (ipip) 7^ - in relativistic field theory this is known as "chiral symmetry 
breaking". Physically the most important outcome is the generation of a gap in the 
fermion spectrum, implying the model describes an insulator. Son postulates that this 



■""In experiments it is only possible to reduce the effective electron charge by mounting the graphene 
layer on a dielectric substrate. 
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insulating phase exists in the corner of the phase diagram corresponding to large e 2 
and small Nf, and in particular that the insulator-conductor phase transition taking 
place at Nf = Nf c in the strong-coupling limit e 2 — > oo is a novel quantum critical 
point. The value of Nf c , and the issue of whether it is greater than of less than the 
physical value Nf = 2, must be settled by a non-perturbative calculation. A recent 
estimate, obtained by a renormalisation group treatment of radiatively-induced four- 
fermion contact interactions, is Nf c = 2.03 |6J. 

The proposed physics is very reminiscent of another 2+lci fermion model, this time 
relativistically covariant, namely the Thirring model with action 



Th 



£ 

a=l 



dxod 2 x 



(5) 



with fi = 0, 1,2, particularly once we insist on units such that vf = 1- Note that in 
contrast to the graphene model the coupling g 2 has mass dimension -1. Once again, 
the model is analytically tractable at large Nf, but exhibits spontaneous chiral sym- 
metry breaking leading to gapped fermions at small Nf and large g 2 [H [S|. Arguably 
the Thirring model is the simplest field theory of fermions requiring a computational 
solution: the location of the phase transition at Nf = Nf c in the strong coupling limit 
has recently been determined by lattice simulation to be Nf c = 6.6(1) [9]. The apparent 
similarity of the two systems has led us to propose a Thirring-like model pertinent to 
graphene, with action 



S 2 = / dx d^x 

a=l 



(6) 



The only difference with (JSJ) is that the contact interaction is now only between the 
time-like components of the fermion current, so that the model is no longer covariant. 

The traditional way to proceed is to introduce an auxiliary boson field V. The 
resulting action 



N f 

S' 2 = / dxod' 

a=l J 



X 



(7) 



reproduces the identical dynamics as ([6]) once V is integrated out. As for (JTJ we assume 
a large- Nf limit to estimate the dominant vacuum polarisation correction; the resultant 
propagator for V is 



l_ Nf |pf 
9 2 



(P 2 



(8) 



In the strong-coupling or large- Nf limits, D 2 coincides with D\ (J3J), implying that the 
fermion interactions are equivalent. It is also the case that linip^oo D 2 (p) = liniA^oo D\(p). 
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This last limit is important because critical behaviour in the Thirring model (J5J) is gov- 
erned by a UV-stable fixed point of the renormalisation group [7j. We aniticipate that 
the model ([6]) is similar and expect its predictions, in particular for critical behaviour 
such as the value of Nf c , to be generally valid for Son's model (CD) in the limit of large A. 

In the following section we present a version of the action ([7]) discretised on a space- 
time lattice, and outline how its dynamics can be investigated by standard Monte Carlo 
simulation techniques. To our knowledge this paper is the first to apply lattice gauge 
theory techniques to graphene. In this first paper we focus exclusively on the equivalent 
of the order parameter (ipip) i n t ne {Nf,g~ 2 ) plane. In Sec. 13.11 we explore the strong 
coupling limit and identify the critical flavor number Nf c , and then attempt to char- 
acterise the transition from insulator to conductor by studying the critical equation of 
state using finite volume scaling. In Sec. 13.21 we switch attention to the physical case 
Nf = 2, and present results from a study of ('i/ji/j) as a function of g 2 . A brief discussion 
of the implications for the graphene model of [5] follows. 



2 Lattice Model and Simulation 

The lattice model studied in this paper is closely related to the lattice Thirring model 
studied in [7j. It is written in terms of staggered lattice fermions, ie. single-component 
Grassmann fields \i X defined on the sites a; of a three-dimensional cubic lattice, by the 
action 



Siatt = ^ Yl XaxV^x [(1 + <5 M oy -jy-e lVx )Xax+fL - (1 + <V>y —e ,v «- 6 )x ax _ / J 

xfia 

XaxXax- (9) 

xa 

The flavor indices a — 1, . . . , N. The sign factors r] Xfl = (— l) a; o+-+x M _i ensure that in the 
long-wavelength limit the the first (antihermitian) term in Si a tt describes the Euclidean 
propagation of Nf = 2N flavors of relativistic fermion described by four-component 
spinors [10J. The fermion mass term proportional to m has been added to provide a 
IR regulator for modes which would otherwise be massless; beyond the usual critical 
slowing down, it is important to stress that simulations directly in the limit m — > 
present severe technical difficulties. The hopping terms in Si a u involve the auxiliary 
boson field V x which is formally defined on the timelike links connecting sites x with 
x + 0. The ^-dependence in the kinetic terms is conventional, and retained to ensure 
continuity with the studies of [HI E]- In order to compare with the formulation of the 
previous section the rescaling g 2 —>■ Ng 2 is required. 

It can be shown that the action (Q can be reexpressed in terms of four-component 
spinors in a form similar but not identical to ([6]). For full details of the relation between 
the two actions we refer the reader to [7j. Here we merely note that the spontaneous 
generation of a condensate (xx) 7^ in the lattice model results in a chiral symmetry 
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breaking pattern U(iV)®U(iV) — >U(N), whereas in the continuum model ([6]) (tpip) 
breaks U(2iVj) — >XJ(Nf)(B)U(Nf) [3]. A term proportional to m explicitly breaks the 
symmetry in either case. It is plausible that the effective global symmetry of the lattice 
model enlarges in the continuum limit, and the correct continuum pattern recovered. In 
what follows we will assume that the chiral symmetry breaking described by (xx) 
is equivalent to the metal-insulator transition. 

The novelty of [9] was the first study of the Thirring model ([5]) by lattice means 
in the strong-coupling limit g 2 — > oo. Since we aim to repeat the strategy here we 
discuss how this was done. First, note that the vacuum polarisation calculation leading 
to the results ( TSTBl) does not go through in quite the same way for the lattice regularised 
model (jUj); rather, there is an additive correction which is momentum independent and 
UV-divergent: 

U latt (p) = Il cont (p) + g 2 J(m), (10) 

where J(m) comes from incomplete cancellation of a lattice tadpole diagram [7]. This 
extra divergence not present in the continuum treatments can be absorbed by a wave- 
function renormalisation of V and a coupling constant renormalisation 



In the large- Nf limit we thus expect to find the strong coupling limit of the lattice model 
at g\ — > oo implying g 2 — > g 2 m = J -1 (m). For g 2 > g 2 im Di att (p) becomes negative, and 
Siatt no longer describes a unitary theory. 

Away from the large- Nf limit, where chiral symmetry may be spontaneously broken, 
there is no analytical criterion for identifying g^ m ; however in this case a numerical 
calculation of (xx) shows a clear peak at g 2 = g 2 cak , whose location is approximately 
independent of both volume and m, indicating an origin at the UV scale [9]. Fig. [I] 
exemplifies this behaviour in the model (Q with Nf = 2 on system volumes L 3 : for 
L > 24 we identify l/g 2 cak — 0.3. Since for orthodox chiral symmetry breaking the 
magnitude of the condensate is expected to increase monotonically with the coupling 
strength, we interpret the peak as the point where unitarity violation sets in, ie. g 2 im ~ 
g 2 eak . In the next section, we shall use simulations performed at g 2 = g 2 eak to explore the 
strong coupling limit, and find evidence for a chiral symmetry restoring phase transition 
at a well-defined Nf c . 

Writing the action in the form XiMijXj, it is possible to integrate out the fermion 
fields analytically to yield the path integral 



Watt 



j VV(detM[V;m}) N . (12) 



Techniques to simulate the physics described by (fl~2l) typically proceed by evolving a 
boson configuration {V} through a fictitious simulation time using a quasi-Hamiltonian 
dynamics in which quantum effects are incorporated via periodic stochastic refreshments. 
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Figure 1: (xx) vs - Vff 2 f° r ^/ = 2. 



We implement this using the hybrid molecular dynamics (HMD) algorithm [TT]. The 
key step in the evolution involves the calculation of a force 



5S_ 

sv 



NtiM' 



_ ± 5M 

w 



(13) 



Since the force can be calculated for arbitrary N, it is possible to simulate the dynamics 
for non-integer N, which is equivalent to regarding the path integral ([121) as the funda- 
mental definition of the model. Of course, only for integer N, and therefore for even 
integer Nf, is it possible to express the theory as a local action in the fermion variables 

In the simulations described in this paper we used a HMD algorithm to perform 
simulations with arbitrary Nf. In principle, this method is not exact in the sense that 
results have a systematic dependence on the discrete timestep Ar used to integrate the 
HMD equation of motion. We have used Ar = 0.0025 on the smallest systems (16 3 , 
24 3 , 32 3 ), Ar = 0.00125 on 16 2 x 48 and Ar = 0.000625 on 16 2 x 64; in all cases we 
checked that the resulting systematic error is smaller than the statistical error. The 
mean trajectory length f = 1.0, and (xx) measured using 10 stochastic estimators after 
every trajectory. Roughly 200-400 trajectories were generated for 16 2 x 48, 64, 600 for 
24 3 ,32 3 and 0(1000) for 16 3 . Further details of the numerical methods used can be 
found in [HE]. 
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3 Results 



We performed simulations on system volumes L 2 S x L t = 16 3 , 24 3 and 32 3 , using fermion 
masses m = 0.01, . . . , 0.04. Because the action (J9j) does not treat spacelike and timelike 
directions equivalently, we also found it useful to explore the consequences of indepen- 
dently varying L s and L t , and thus in addition studied 16 2 x 48,64 ; 24 2 x 32,48; and 
32 2 x 24. As we shall see, the aniostropic nature of the model's dynamics results in 
the systematic effects due to finite L t being much more important than those due to 
finite L s . The only observable discussed in this initial study is the chiral condensate 
(XX) = (trM- 1 ). 



3.1 Strong Coupling Limit 
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Figure 2: l/^ cak vs. N f . 



As described above, to explore the strong coupling limit g\ —>■ oo we made the ansatz 
g 2 = (?p Cak , where g 2 eak denotes the location of the peak in (xx) a ^ given Nf. Fig. [2]shows 
Vflpeak(-^/) ^ or some representative lattice volumes and fermion masses, confirming that 
its value, arising as it does from UV lattice artifacts, is to good approximation volume 
and mass independent. The behaviour is qualitatively similar to that found for the 
strong-coupling Thirring model shown in Fig. 3 of [9j. We see that l/(?p eak decreases as 
Nf increases from 2 to 4.75, at which point the curve reaches a minimum. There is then 
a steep increase at Nf ^4.9 followed by a levelling off, implying a significant change 



7 



in the model's strong coupling behaviour. Our interpretation, to be supported below 
by a study of the equation of state, is that for Nf below the change the model is in a 
chirally broken phase, and that above the change chiral symmetry is restored, implying 
lim m ^o(xx) = 0- Using this criterion we identify the critical flavor number for chiral 
symmetry restoration in the strong coupling limit, with a conservative error, as 



N 



fc 



4.8(2). 



(14) 



0.25 



0.2 



0.15 



0.1 



0.05 



u 5 



° 24 x24 
24 2 x32 
• 24 2 x48 
D 16 2 x48 



_i_ 



N f 

Figure 3: (xx) vs - Nf at m — 0.01 on various lattice volumes. 

In the rest of this section we analyse data taken at g 2 = g^ eak in an attempt to 
determine the critical equation of state (xx( m > Nf)) m the strong coupling limit. Note 
that typically 4-6 independent simulations were used to identify <?p eak for each Nf. In 
the thermodynamic zero temperature limit, a simple ansatz for the approximate scaling 
behaviour close to a quantum critical point is given by 



m = A(N f -N fc )(xx)* + B(xx) i 



(15) 



whence (xx)\ 



N f =N fc 



(x m* and the conventional exponent in (xx)\m=o oc {Nf c — Nf) 13 
for Nf < Nf c is given by (3 = (5 — p)^ 1 . We must, however, take account of the fact 
that our data is taken on finite systems. Fig. [3] shows data for m = 0.01 from various 
lattices: a comparison between 16 2 x 48 and 24 2 x 48 demonstrates that the dominant 
finite volume effects are due to varying L t , while the effects of finite L s are negligible for 
L, > 16. 
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The theory of finite volume effects in models such as (Q with anisotropic correlations 
is outlined in [12]. Near a critical point it is possible in principle to distinguish two 
correlation lengths £ s and £t, each diverging with a distinct critical exponent v s , v t 
as Nf — * Nf c . In d spacetime dimensions these are related to conventionally-defined 
exponents governing scaling of the order parameter and its associated susceptibility via 
a generalised hyperscaling relation 

u t + (d-l)v s = 7 + 2/5. (16) 

Motivated by Fig. [3J in our analysis we take a pragmatic approach and assume all volume 
effects are due to finite L t . The ansatz for the modified equation of state, inspired by a 
renormalisation group analysis [7] is then 

m = A[(Nf - N fc ) + CL t \xxY + B(xx) S ; (17) 

we fit this form to our dataset with a least squares fit. Our complete dataset contains 124 
data taken at various Nf, m, L s , L t (recall the value of l/<7p eak must be independently 
determined for each parameter set, so the simulation effort involved is considerable; 
approximately 100 000 processor hours using 2.4GHz Opterons were required). 

Experience with previous models shows that the fitted equation of state is very 
sensitive to assumptions made about the scaling window (ie. the ranges of Nf and m 
to include in the fit), and the smallest volume to include in the scaling ansatz (ITTj) . For 
this reason we judge it best to present a compilation of different fits in Table [TJ We 



fit 


# 


A 


B 




S 


P 


C 




xVdof 


Power, 


28 


0.31(3) 


41.5(15) 


3.81(3) 


3.96(3) 


0.87(3) 






4.8 


Power, m > 0.02 


18 


0.30(7) 


87(55) 


3.87(9) 


4.4(4) 


0.82(6) 






5.9 


Power, m > 0.03 


12 


2.1(10) 


3800(18) 


4.3(1) 


6.0(1) 


1.3(1) 






6.4 


FVS, m = 0.01 


53 


1.5(7) 


63(22) 


4.60(15) 


3.9(1) 


1.3(1) 


9.7(10) 


1.7(2) 


4.0 


as above L t > 24 


48 


4.5(32) 


161(97) 


4.95(16) 


4.0(1) 


1.6(2) 


7.9(5) 


2.1(2) 


4.0 


as above Nf < 4.5 


46 


712(100) 


2.3(3)xl0 4 


4.46(9) 


5.25(4) 


3.00(4) 


15(3) 


1.2(1) 


5.4 


FVS, all m 


96 


0.23(1) 


19(2) 


3.85(4) 


4.03(8) 


0.88(1) 


17.5(17) 


1.10(5) 


6.3 


as above L t > 32 


70 


0.20(1) 


10.5(1) 


4.7(3) 


3.6(1) 


0.82(2) 


6.0(6) 


2.6(8) 


5.1 


as above Nf > 3 


60 


0.21(1) 


237(106) 


4.6(7) 


5.5(3) 


0.86(2) 


8.1(34) 


2.1(1.1) 


3.1 


N f > 3 L t > 24 


75 


0.21(1) 


352(137) 


4.1(1) 


5.8(2) 


0.87(2) 


15(4) 


1.3(2) 


3.4 


FVS, m > 0.02 


43 


0.19(2) 


10.4(18) 


3.76(11) 


3.55(13) 


0.78(3) 


24(15) 


0.9(2) 


7.3 


as above L t > 32 


32 


0.16(2) 


6.5(14) 


3.9(4) 


3.21(15) 


0.72(4) 


12(19) 


1.2(9) 


6.4 



Table 1: Various fits to the Equations of State (H]) "Power" and (TT7]) "FVS 



tried fits to both the "Power" equation of state ([15]) . with 5 free parameters, using data 
from a single lattice size 16 2 x 64, and fits using the finite-Lj "FVS" scaling form (ITTj) 
with 7 free parameters^ In the latter case data with Nf > 5 was excluded from the fit 

2 Note it is not possible to use hyperscaling to constrain the value of ft, as done in [71 [HI [3] • 
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because their small errorbars destabilised the fits; since this data probably comes from 
the chirally symmetric phase there may be a small systematic error in the identification 
of (?p Cak across the transition. 

Fits to (fT5i) favour Nf c ~ 3.8 — 4, and p ~ 0.9. These values are also favoured by the 
most comprehensive FVS fit to the 96 datapoints with Nf < 5. There is no evidence that 
discarding m = 0.01 data, which may be most prone to finite volume artifacts, improves 
any of the fits. On the other hand, discarding L t = 16 and perhaps L t = 24 does have 
a significant effect on the fitted values of Nf c , p and v t in the FVS fits. In these cases 
the fitted 5 ~ 4. However, once data with extremal values of Nf is excluded, on the 
assumption that they lie outside the scaling window, the fitted values of S rise to ^ 5. 
In almost all cases the fitted value of v t exceeds 1, though often not by a statistically 
significant margin. 

16 2 xl6ra=0.01 
O 24 2 x24 ra=0.01 

• 32 2 x24 m=0.01 
24 2 x32 m=0.01 
32 2 x32 m=0.01 

° 16 2 x48m=0.01 

• 24 2 x48 m=0.01 

• 16 2 x64m=0.01 
n 24 2 x24 m=0.02 

32 2 x32 m=0.02 
n 16 2 x48m=0.02 
■ 16 2 x64ra=0.02 
O 24 2 x24 m=0.03 
« 32 2 x32 m=0.03 
O 16 2 x48m=0.03 

• 16 2 x64m=0.03 
A 16 2 x64m=0.04 



"3 4 5 6 7 8 

N f =N f +C/L~ 1/w < 

Figure 4: Finite volume scaling fit to (fT7|) to data with m — 0.01 (circles), 0.02 (squares), 0.03 (dia- 
monds) and 0.04 (triangles), in terms of NL 

Our favourite fit, yielding the smallest x 2 /dof, emerges from the 60 datapoints with 
Nf G [3,5) and L t > 32. Another reason for preferring this is that the fitted Nf c is 
consistent with the value (JT3J) coming from the behaviour of <7p Ba k(-W/)> which could be 

regarded as an additional constraint on the global fit. The fit is plotted in Fig. H] in 

i_ 

terms of the control parameter in the thermodynamic limit N'f = Nf + CL t , so that 
data with differing L t should collapse onto a single curve for each value of m. 

To summarise: this "best" fit provides a reasonable description of the data in the 
window 4.5 ^ N' f ^ 6, in particular for the smallest mass m = 0.01; fits of the form 
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n 



0.2 



o.i 



ft. 
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([IT]) are capable of yielding a fitted Nf c consistent with (TH|) ; and the preferred value of 
Ut « 2 once this consistency criterion is applied. 



3.2 Nf = 2 

Next we turn our attention to the physical case iV/ = 2. Since Nf c > 2, we expect chiral 
symmetry to be broken at strong coupling and potentially restored at some finite g 2 R . 
Accordingly we study (xx) as a function of l/g 2 and here study additional values of the 
mass parameter m. Our results are summarised in Fig [5j 
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Figure 5: (xx) vs - Vff 2 f° r Nf = 2. 

As before, we have attempted to fit a critical equation of state using the forms ( fl5l) 
to data from 16 2 x 64, and ( IPTI) . in each case replacing (Nf — Nf c ) with (g~ 2 — g~ 2 ). No 
stable fits were found unless data with l/g 2 1 were excluded. Our most comprehensive 
fit, using the FVS form and restricting the data to L t > 48, l/g 2 < 0.9, yields l/g 2 = 
0.632(6) and is shown in the figure. All the fits we found identify a l/g 2 ~ 0.6 ^> l/#p eak , 
but all clearly undershoot the data at weaker couplings by a considerable margin. We 
conclude that the eqution of state ansatze (I15|I17I) are inadequate to describe the data. 

Instead, we will distinguish between a "strong coupling" regime l/g 2 ~ 0.75 where 
(xx) is numerically large and finite volume effects are negligible and a "weak coupling" 
regime l/g 2 ~ 0.75 where the opposite holds true. Two comments about the weak cou- 
pling regime are worth making. First, as is clear from Figs. [5] and El finite volume effects 
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are unexpectedly large, and indeed increase in relative importance until l/g 2 ^ 1; it is 
this feature which has made a global FVS fit impossible. In a conventional symmetry- 
breaking scenario by contrast one expects the finite volume effects to be larger in the 
broken phase, where there are long-range correlations due to Goldstone bosons. Sec- 
ondly, in a chirally symmetric phase one expects (xx) <X- m f° r small m and weak in- 
teractions; inspection of the data plotted in Fig. [6] appears to imply either that a linear 
extrapolation to m — > would yield a non-vanishing order parameter, or alternatively, 
that chiral symmetry restoration for these values of g 2 < g 2 would require (xx{ m )) f° 
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Figure 6: (XX) vs - 771 f° r Nf = 2 in the weak-coupling regime. 

exhibit some negative curvature, for which there is some tentative evidence in the figure. 
We conclude that either chiral symmetry remains broken at weak coupling, or that there 
are long-range correlations in this region. 

We briefly consider alternative scaling scenarios. In the chiral and thermodynamic 
limits two other kinds of behaviour are possible to envisage, and are not currently ex- 
cluded by our data. Firstly, a form (xx) = Ae~ B / 9 would predict broken chiral symme- 
try for all g 2 . This is barely credible, since the x - X forces in this model are weaker than 
for the Nf = 2 Thirring model, where two independent simulation studies find a second 
order chiral restoring transition at l/g 2 — 1.9 [3 |T3j (Fig. [5] also shows Thirring data 
from 16 3 [7j). Second, (xx) = Ae~ B ^ 9c ~ 9 2 ^ q describes chiral symmetry restoration 
via an infinite order phase transition at g 2 = g 2 . Without a reliable finite volume scal- 
ing hypothesis we cannot estimate g 2 , B or q. By analogy with the Kosterlitz-Thouless 
transition in 2d systems, though, this scenario predicts a critical and hence strongly fluc- 
tuating system for g 2 > g 2 which could plausibly account for the observations reported 
above. 
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4 Discussion 



In this paper we studied a model ([6]) which has very similar properties, including the 
same global symmetries, as the graphene-related model rececently proposed to describe 
quantum critical behaviour in the (Nf, e 2 ) plane [5]. Using a simulation strategy devised 
for the 2+ld Thirring model [9], we have identified the critical number of flavors separat- 
ing insulating from conducting phases in the strong coupling limit to be Nf c = 4.8(2). 
This implies that the strong coupling limit of graphene with Nf = 2 is an insulator. 
Since the properties of a critical point should be universal, we expect this result to be 
a robust prediction of our work, thus furnishing the first systematic non-perturbative 
prediction of quantum critical behaviour in this system. We also managed a reasonable 
fit of our strong-coupling data to an equation of state describing a continuous phase 
transition at the critical point, and obtained estimates for the critical exponents. 

The fitted value of v t gives partial information about the nature of correlations in the 
vicinity of the fixed point. Substituting our favoured values 5 = 5.5, p = 0.86, u t = 2.1 
into (fT6l) . we obtain v s = ^7 — 0.83. Without studying the order parameter susceptibility 
we have no independent estimate of 7, but note that it would need to have a value of 
0(6) in order for v s to exceed v t . Now, there is a further relation governing the scaling 
of critical correlation functions [12] : 



where order parameter correlations (xx(fyxx( x s,t)) oc x s J a,t at criticality, with the expo- 
nent taking the appropriate value depending on whether the displacement x is timelike 
or spacelike. The ratio r] s /i] t > 1 if v s /v t < 1 and vice verso,. However, Tj s /Tji may 
be identified with the dynamical critical exponent z characterising the quantum critical 
point, in the sense that the dynamics remains invariant under the scale transformation 
x — > Ix; Xq —>■ £ z Xq. By considering the anomalous dimension of the Fermi velocity 
using the 1/Nf expansion, Son [5] has obtained z < 1, which has implications for the 
stability of the quasiparticle excitations. If we assume that the same critical exponent 
governs both quasiparticle and order parameter correlations, then reconciling the two 
calculations requires an unusually large value of 7. 

Next, we set Nf to the physical value 2 and studied the chiral order parameter as 
a function of coupling strength. Here our results are harder to interpret; we observe 
a crossover from strong- to weak-coupling behaviour at 1/g 2 — 0.75, but were unable 
to model the equation of state, leaving the nature of the weak-coupling regime un- 
clear. There is evidence, both from the large finite volume effects and the curvature in 
(xx(m)), for strong correlations. Work is currently in progress to study the quasiparticle 
propagator in order to explore the dispersion relation and expose any quantum critical 
behaviour from an independent direction, and also to further investigate the nature of 
the fluctuations in the weak-coupling phase at Nf = 2. 



{d-2 + r] t )v t = (d-2 + r] s )u s => 



vt_ = l + r) s 
v s 1 + Vt ' 
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